#!/usr/bin/python

import math
from math import log
from math import pow
from math import pi

from os import system

c=2.998e8# m/s
########################################################################
###### Parameters ######################################################
########################################################################
distance = 30000 # m

points_nsc = [(0,0,0),
			(0,0,1),
			(0,0.7,1.7),
			(0,0.7,2.7)]
			
points_units = 2e3 # m
# Applying points scale
points = list()
L=0
lengths = list()
for i in range(0,len(points_nsc)):
	coords = points_nsc[i]
	points.append((coords[0]*points_units, coords[1]*points_units, coords[2]*points_units))
	if i > 0:
		L += pow(pow(points[i][0] - points[i-1][0],2)+pow(points[i][1] - points[i-1][1],2)+pow(points[i][2] - points[i-1][2],2), 0.5)
	lengths.append(L) 

# L - summary length
# lengths[node_i] - Length of all segments previous to this node


gamma = -log(0.7)/(L/c) # Approximately 0.5 speed decrease [Hz]
v0 = 0.7 # [light speed]

def delay(l):
	tau = -1/gamma*log(1-l*gamma/(v0*c))
	return tau

def v_c(l):
	return v0 - gamma*l/c

########################################################################
########################################################################
########################################################################
names = list()
maxksi = 4
for ksi in range(0,maxksi+1):
	for field in ["E", "H"]:
		alpha = ksi/float(maxksi) * pi/2*0.9
		if field == "E":
			n = [-math.sin(alpha), 0, math.cos(alpha)]
			yunit = "E[V/m]"
		else:
			n = [0,1,0]
			yunit = "H[A/m]"
		obs = [distance*math.cos(alpha), distance*math.sin(alpha)]
		deg = round(alpha/pi*180)
		name = "alpha %i deg"%(deg)
		filename = "alpha_%ideg%s"%(deg, field)
		names.append((field, name, filename));
		
		xml = open("in.xml", "w")
		xml.write("<model>\n")
		xml.write("	<source>\n")
		for i in range(0,len(points)-1):
			t_start = delay(lengths[i])
			p0 = points[i]
			p1 = points[i+1]
			xml.write("		<discharge t_start='%e'>\n"%(t_start))
			xml.write("			<current>\n")
			xml.write("				<shape type='expexp' alpha='1e4' beta='1e5' Q='3.87'/>\n")
			xml.write("				<delay type='exponential' v_c='%e' gamma='%e'/>\n"%(v_c(lengths[i]), gamma)) # Doing so is OK, because \frac{dv}{dt} = -gamma*v - matching V is enough to reproduce the process
			#xml.write("			<decrease type='exponential' eta='0.7e-3'/>\n")
			xml.write("				<decrease type='none'/>\n")
			xml.write("			</current>\n")
			xml.write("			<trajectory type='line' r0='[%e,%e,%e]' r1='[%e,%e,%e]'/>\n"%(p0[0],p0[1],p0[2],p1[0],p1[1],p1[2]))
			xml.write("		</discharge>\n")
		xml.write("	<mirror/>\n")
		xml.write("	</source>\n")
		xml.write("	<observer r='[%e,0,%e]' field='%s' n='[%e,%e,%e]'/>\n"%(obs[0],obs[1],field, n[0], n[1], n[2]))
		xml.write("</model>\n")
		xml.close()
		system("cat in.xml | ./qgr > %s.dat; cp in.xml %s.xml"%(filename, filename))


cfg=open("cfg.m", "w");
cfg.write("signals = {")
first = True
for row in names:
	if first:
		first = False
	else:
		cfg.write(";\n")
	cfg.write("'%s', '%s', '%s'"%(row[0], row[1], row[2]))
cfg.write("};\n")
cfg.close()
system("octave vis.m")
print "Done"


